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Abstract 

An electronically coarse-grained model for water reveals a persistent vestige of the liquid-gas transition 
deep into the supercritical region. A crossover in the density dependence of the molecular dipole arises 
from the onset of non-percolating hydrogen bonds. The crossover points coincide with the Widom line in 
the scaling region but extend further, tracking the heat capacity maxima, offering evidence for liquid- and 
gas-like state points in a one-phase fluid. The effect is present even in dipole-limit models suggesting that 
it is common for all molecular liquids exhibiting dipole enhancement in the liquid phase. 
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There is a distinct boundary between liquid and vapor phases of a substance separated by a 
coexistence line over a finite range of pressures (p) and temperatures (T). Crossing the line results 
in a discontinuous density change - the hallmark of a first-order phase transition. The coexistence 
line terminates at a critical point in the (T, p) plane beyond which the thermodynamic distinction 
between liquid and gas phases is lost and a single-phase supercritical fluid is formed. In this 
region, there are no further phase boundaries and the supercritical fluid permits a continuous path 
from liquid to gas over a broad range of thermodynamic states [1]. In water, for example, this 
tunability can be exploited to control solvation properties making supercritical water an important 
ingredient in industrial processes and a promising “green” alternative to chemical solvents [2-4]. 

While the classical thermodynamic picture is well understood, the molecular nature of super¬ 
critical fluids remains a largely open question. Key issues for water are (1) structure and hydrogen 
bonding (HB) in the supercritical region; (2) the electronic redistribution which occurs in the water 
molecules as the HB network is reconfigured; and (3) most broadly, the relationship of the super¬ 
critical fluid to fundamental liquid and gaseous states of matter that exist as separate phases below 
T 

In this context, the observation that certain thermophysical response functions exhibit max¬ 
ima in the transition region between gas-like and liquid-like states of the supercritical phase has 
prompted the suggestion that a remnant of coexistence may extend into the supercritical region 
[5]. The locus of such maxima is referred to as the Widom line [6]. Direct confirmation of dis¬ 
tinct gas-like and liquid-like state points, however, is lacking though some evidence of dynamical 
crossover phenomena has been reported in noble gases [5, 7]. Analogous concepts may apply at 
low temperatures where a Widom line is proposed to extend from a hidden liquid-liquid critical 
point [6] influencing the properties of supercooled water. 

Supercritical water poses unique challenges for molecular models because conventional em¬ 
pirical interaction potentials capture neither the electronic redistribution which occurs along an 
isotherm within the supercritical region of the phase diagram (dipole, quadrupole and higher 
manybody polarization responses) nor the manybody dispersion forces which have emerged as 
an important aspect of water physics [8, 9]. Moreover, as these empirical models are typically op¬ 
timized to describe liquid water near ambient conditions, their transferability to the supercritical 
fluid is questionable and the resulting physical insight into the molecular processes which occur 
above the critical point is limited. 

A new approach for materials simulation has been proposed recently wherein the molecular 
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electronic distribution is not represented by a parametrized mean-field model, but rather by a 
coarse-grained description based on embedded harmonic oscillators - the quantum Drude oscil¬ 
lator (QDO) model [10, 11]. As the oscillators are treated quantum mechanically, this model 
generates the complete hierarchy of many-body polarisation and dispersion phenomena which 
when solved in strong coupling, as here, ensures that all potential symmetry-breaking interactions 
are present within Gaussian statistics. The responses (to all orders) include, but are not limited 
to, inductive responses leading to distortions of the electronic charge distribution and to van der 
Waals and higher-order dispersion interactions arising from quantum mechanical charge density 
fluctuations. The model parameters are fit to leading-order responses of the isolated molecule 
providing a perfect low-density limit [9, 12, 13]. Here we apply the model with the short-range 
repulsion parameterized to a high-level ab initio calculations of the water dimer, (see Refs [9, 14] 
for model details). 

In prior work [9], using the methodology described in [10-13] to simulate the model, we have 
determined the location of its critical point as (T c = 649(2) K, p c = 0.317(5) g/cm 3 }, which is 
in accord with experiment (T c = 647.096 K, p c = 0.322 g/cm 3 } [15]. The dielectric properties 
along coexistence as well as the orthobaric densities of both the gas and liquid branches are well 
described [9] . We regard the prediction of an accurate critical point as a prerequisite for reliable 
examination of the supercritical phase. We are aware of no other description able to predict an 
accurate critical point from the properties of the isolated molecule. 

Having established the critical point and the dielectric properties of the model, we are now in 
position to examine the molecular properties across several supercritical isotherms and to draw 
broader conclusions about the molecular nature of supercritical water. To set the context for this 
part of the work, we draw attention to recent studies which have explored the extent to which 
the notion of liquid-gas coexistence can be extended past the critical point. Typically, this ex¬ 
trapolation is based on identifying the loci of maxima in second-order thermodynamic response 
functions such as heat capacity (C p ), thermal expansivity (a) or isothermal compressibility (k t ) 
along isothermal paths. In atomic systems (noble gases), evidence of dynamic anomalies near the 
points of maxima in C p has been reported [7] and interpreted as a signature of liquid- and gas-like 
behavior. However, it is clear that the extrema of different thermodynamic functions rapidly di¬ 
verge from each other above the critical point and only collapse onto a universal Widom line in 
the asymptotic scaling region near the critical point as described by a scaling theory developed in 
Ref. [1]. Beyond this region, and further into the one-phase supercritical field, it is unclear whether 
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any specific thermodynamic observable separates distinct regimes which can be unambiguously 
associated with liquid- or gas-like behavior. 

We now consider the dipole moment of water as a local reporter of molecular environment. In¬ 
dividual water molecules which are adaptive and responsive to their environment, become strongly 
polarized as a result of the highly directional bonding of the liquid-like state [16]. We explored the 
supercritical region along several isotherms from 673 K to 1083 K and for densities from 2 mol/l 
to 60 mol/l spanning from 10 MPa to 2.4 GPa. The results of our calculations of the molecular 
dipole at two isotherms in the supercritical region of the phase diagram are shown in Fig. 1. The 
dipole’s density dependence exhibits a change in slope at p = 17.1 mol/l at 673 K and 18.9 mol/l 
at 873 K. The corresponding densities of heat capacity maxima, obtained from the IAPWS-95 
reference equation of state for water [15], are 17.25 mol/l and 18.91 mol/l. The dipole crossover is 
still detectable at T = 1083 K = 1.7 x T c above which it disappears in statistical noise. It becomes 
more pronounced with decreasing temperature and below T c it occurs in the two-phase region. 

In order to investigate further the crossover in the dipole moment we analyze another micro¬ 
scopic discriminator between the liquid-like and gas-like states - hydrogen bonding (HB). Al¬ 
though the presence of a HB network in the supercritical region is well established by various 
techniques [17-19], the structure and the topology of the network in liquid water is a matter of 
discussion [19-22]. Examining HB connectivity, we find that local tetrahedral order persists down 
to low gas-like densities at supercritical temperatures. Using the distance-angle geometric defini¬ 
tion of the hydrogen bond of Ref. [23] we observe a similar, albeit less pronounced, crossover in 
density dependence of the average number of hydrogen bonds per molecule, which occurs at the 
same densities as the crossover in dipole moment, shown in Fig. 1. 

The number of hydrogen bonds per molecule at 673 K and density p = 30 mol/l (0.54 g/cm 3 ) 
is 40% of that of ambient water, 300 K, 55.32 mol/l, where n H B = 3.71. This can be compared 
with 29% of ambient at the same temperature and 28.9 mol/l estimated from NMR chemical 
shift measurements [17]. The number of hydrogen bonds per molecule at the crossover density, 
which decreases slightly from 1.08 at T = 673 K to 0.85 at T = 1083 K, is well below the 
percolation threshold which occurs at n H B = 1.53(5) according to Ref. [24]. This observation is 
in accord with Raman experiments [22], which are interpreted to indicate no tetrahedral hydrogen 
bonding at the critical point. We note that collective properties which may show dynamic or 
viscoelastic anomalies may also be detectable by x-ray scattering techniques or other experimental 
probes of relaxation processes. We therefore conclude that the onset of the formation of minimal 
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HB associations corresponds to the dipole crossover point in the supercritical field which, with 
decreasing temperatures, evolves to the liquid-gas coexistence line in the two phase region. 

The results of the study are summarized in Fig. 2, where we plot the positions of the observed 
crossover points in the (T, p) plane along with the portion of the coexistence line ending in the 
critical point and experimental data showing maxima in various response functions. Shown here 
are specific heat, C p , isothermal compressibility, k t , and thermal expansivity, a. Asymptotically 
close to the critical point, in the region T c < T < T w , the dipole crossover points occur along 
the Widom line, where the loci of all response function maxima coincide. Moreover, the dipole 
values are unambiguously vapor-like on the low-p side and liquid-like on the high -p side. On the 
molecular scale it therefore appears possible to partially extend the notion of coexistence and the 
liquid-gas transition well into the thermodynamically single-phase fluid as a remnant. 

At higher temperatures, in the non-asymptotic region T w < T < T H , the loci of response 
function maxima diverge from the Widom line but it appears that the dipole crossover most closely 
tracks C P (T ). Above T H the dipole derivative transition becomes undetectable and fluid appears 
homogeneous at all lenthscales. We note that the formation of HB dimers creates an additional 
thermal reservoir which contributes to the specific heat in liquid water. We therefore speculate 
that this may be the molecular origin of the close link revealed here between dipole-crossover, 
hydrogen bond formation, remnant liquid-gas coexistence behavior in the supercritical phase and 
the locus of maxima in C p in and beyond the scaling region. No such direct link exists between 
HB formation and thermal expansivity or isothermal compressibility the maxima of which follow 
very different trajectories out of the scaling region. 

Finally, we examine the structure and density at several thermodynamic state points for which 
experimental data are available [25]. We focus first on the T = 673 K isotherm and compare the 
results to neutron diffraction measurements [25]. Predicted densities along this isotherm are shown 
in Fig. 3 together with data obtained from IAPWS-95 reference equation of state [15]. Since our 
QDO model does not involve any condensed-phase parameterisation, the agreement arises as a 
prediction and gives confidence in the previously presented results. 

Further structural data is presented in the insets to Fig. 3, which show three-dimensional first 
coordination shell plots for oxygen (red) and hydrogen (white) of surrounding water molecules. 
The isosurfaces are drawn at 2x (bulk density) levels. These illustrate the extent of hydrogen bond 
reconfiguration which occurs at various densities along the isotherm. The partial radial distribution 
functions are shown in Fig. 4 for T = 673 K, p = 340 MPa (red lines) in comparison to the neutron 
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diffraction data of Ref. [25] (blue lines). Subtle features of the experimental H-H distributions are 
present in the model, and the 0-0 correlations are captured as well. 

In Fig. 1, the mean molecular dipole moment \p,\ is seen to vary linearly with density in two 
regions, which can be understood using the following argument: The induced dipole moment is 
proportional to the local field acting on a water molecule, which can be decomposed into contri¬ 
butions from the first coordination shell and from the rest of the fluid. As ab initio and induction 
model calculations of small water clusters [26] have shown, the water dipole moment scales lin¬ 
early with the cluster size for up to five members in the cluster. A mean field estimate assuming 
the average distance to scale like p ~ 1 /3 leads to the same conclusion. More deeply, the Onsager 
reaction field strength scales linearly with density [27]. 

For strongly associated liquids like water, the Widom line has been connected with a percolation 
transition in the hydrogen-bond network [28]. However, our results support the view that the 
origin of the crossover and the extrema in the response functions is more local. We show that 
this behavior extends much further into the notionally one-phase region and can be detected as 
far as T = 1.7 x T c tracking most closely the maxima in the specific heat. Only along the 
isotherms beyond these temperatures does it appear that supercritical water is a homogeneous fluid 
on the molecular scale exhibiting no detectable remnants of the liquid-gas transition. The extent 
of hydrogen bonding across the Widom line and beyond the scaling region can, in principle, be 
indirectly accessed by Raman scattering measurements of the OH stretch band which is a reporter 
of hydrogen bonding and the integrity of the network. Collective properties which may show 
dynamic or viscoelastic anomalies may also be detectable by x-ray scattering techniques or other 
experimental probes of relaxation processes. 

In order to address the question of generality of the phenomena reported here we studied two 
dipole polarizable models (polarization treated in the electric dipole limit), namely the classi¬ 
cal Drude oscillator model, SWM4-NDP [29], for which the critical parameters are known [30]: 
T c dip = 0.89T C , p dip = 0.90p c , p dip = 1.02p c , in terms of experimental values for water. The second 
is a polarizable version of the Stockmayer model, the ‘minimalist’ (purely dipolar) fluid model 
with a different local association topology [31], and estimated critical parameters T c s = 1.1T C , 
p s c = 1.5 p c , p s c = 1.02p c (see Ref. [14]). Finite size effects are investigated and ruled out in 
Supplementary Material [14]. The results for three models, presented in Fig. 5, show that the 
dipole moment change is most pronounced in fully responsive QDO, but also emerges from dipole- 
polarizable water model and the polarizable Stockmayer model. The reduced crossover density for 
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QDO model p* = 0.90, and shifts to 0.63 for SWM4-NDP and to 0.56 for the Stockmayer model. 
Therefore, our study has revealed an intriguing new property of associating liquids, the extension 
of gas-liquid critical effects on the molecular dipole moment. 

In summary, we find direct evidence of molecular-scale heterogeneity in supercritical polar 
fluids. A clearly identifiable transition between dissociated (gas like) and associated (liquid like) 
regimes is evidenced by a change in the density dependence of the molecular dipole moment ac¬ 
companied by a transition in the hydrogen bond connectivity. This characteristic signature has 
been observed in a fully responsive electronic model of water which provides an excellent predic¬ 
tion of the ambient, critical and supercritical properties but is also present in simpler dipole-limit 
models with two different association topologies - the minimal model class in which the the be¬ 
havior may be observed. While the simpler models do not predict the critical point accurately for 
water they nevertheless provide compelling evidence that the phenomenon is likely to be general 
for polar liquids. The observations provide a molecular basis for the variety of response function 
anomalies which define the “Widom line” asymptotically close to the critical point. This is also 
significant outside the scaling regime where the maxima of thermodynamic response functions 
follow different trajectories making it unclear which, if any, separate liquid-like from gas-like re¬ 
gions. Here we observe that the dipole crossover points follow only the heat capacity maximum 
outside of the scaling region and their correlation to rudimentary hydrogen bonding establishes an 
unambiguous extension of the liquid-gas coexistence line deep into the supercritical phase. 

This work was supported by the NPL Strategic Research programme, the Engineering and 
Physical Sciences Research Council (EPSRC) and the European Metrology Research Programme 
(EMRP). Generous allocation of time on BG/Q at STFC Hartee Centre, UK, is gratefully acknowl¬ 
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FIG. 1. (Color online). Average molecular dipole moment as a function of density at T = 673 K, blue 
diamonds, and at T = 873 K, orange circles. The triangle denotes the isolated monomer value (parameter 
of the model). Green squares and right scale, average number of hydrogen bonds per molecule as a function 
of density at T = 673 K. The lines are the linear fits to data. Insets illustrate the variation in electron density 
with pressure, where the pink and blue isosurfaces correspond to gain and loss of the electron density. 
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FIG. 2. (Color online). (T, p) phase diagram of water in the proximity of the critical point in units of critical 
temperature, T c , and critical pressure, p c . The liquid-vapor coexistence curve (pink), the loci of maxima in 
heat capacity (green), thermal expansivity (orange), and isothermal compressibility (blue), obtained using 
the IAPWS-95 reference EOS [15]. The blue diamonds represent our results derived from the analysis of 
Fig. 1. The black dot marks the critical point. The inset, the corresponding (p, T) phase diagram in reduced 
units, shows the gas and liquid branches of the coexistence curve, the locus of the heat capacity maxima, 
and superimposed are linear fits to the dipole moment (the blue dashed lines) in gas- and liquid-like regions 
with vertical offsets to map the corresponding temperatures at the crossover. The crossover densities are 
denoted by the diamonds. The scaling theory of Ref. [1] predicts the “true Widom line” for water where 
dp/dT > 0 follows kt- 
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FIG. 3. (Color online). Density of water along the T = 673 K, 773 K, and 873 K isotherms (from top to 
bottom) as obtained in APIMD-QDO simulation (symbols) vs experimental data(full lines). All simulated 
states, apart from the first four in each isotherm, are in the conventional supercritical region [2]. The dashed 
line indicates the crossover densities. Insets, the 3D plots showing first oxygen (red) and hydrogen (white) 
coordination shells, illustrate remarkable persistence of tetrahedral structure at T = 673 K. 
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FIG. 4. (Color online). The partial radial distribution functions obtained in the QDO water model simulation 
under supercritical conditions at T = 673 K and p = 340 MPa (red lines). Also shown for comparison are 
distribution functions obtained from neutron scattering [25] at the same thermodynamic conditions (blue 
lines). 
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FIG. 5. (Color online). Relative (to the gas-phase value) enhancement of the molecular dipole moment 
as a function of the reduced density. Red diamonds, quantum Drude model; green circles, classical Drude 
oscillator model (dipole limit); blue triangles, polarizable Stockmayer model (dipole limit). All calculations 
were performed at a reduced temperature T* = 1.036 T c in terms of corresponding critical temperature T c . 
Straight lines are to guide the eye. 
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The layout 

A description of the model and simulation methods are provided first. Second, the hydrogen bond 
definition employed in the text is described and justified. A mean-field treatment of the scaling of molecular 
moment with density is given next. The Stockmayer model is described next followed by a discussion 
of finite size effects. Last, we address the different local association topology of water models and the 
Stockmayer fluid. 


The model and simulation details 

For water, which has nearly isotropic polarizability, we follow the model construction reported in 
Ref. [1]. The quantum Drude oscillator (QDO) model for water, sketched in Fig. 1, is based on rigid 
TIP4P geometry of the monomer (roH = 95.72 pm, ZF10F1 = 104.52°). A drudon (quantum Drude os¬ 
cillator) is tethered to a point M along the ZHOH bisector at a distance tom = 26.67 pm from the oxygen 
nucleus. The model parameters are chosen such that the long-range responses match reference values of the 
monomer and dimer (gas-phase polarizabilities and pair-dispersion coefficients). Three fixed point charges, 
two on hydrogens, q\\ = 0.605|e|, and one on M-point, q\i = —2 qu, are taken to represent the isolated 
molecule limit, matched to the low order electrostatic moments. The drudon is also tethered to M-point. Its 
parameters in atomic units are: mass, mo = 0.3656 m e , charge, q\) = — 1.1973|e|, and angular frequency, 
ui = 0.6287 Eh/ft. 

The Coulomb forces are regularized at short range by assigning small Gaussian widths to the charges, 
<th = (j£) = O.lao, and gq = 1.2ao for the QDO center. To correct for missing in the QDO model Fermion 
exchange a biexponential term 


4>{r ) = Ati exp(—Air) + k- 2 exp(—A 2 r) 


has been added between oxygens, providing the required repulsion. The parameters k = {613.3,10.5693} 


2 


and A = {2.3244,1.5145} were obtained from a fit to the difference between a reference quantum chemical 
dimer ground state potential energy surface (CCSD(T)-level) and that of the QDO’s, calculated using norm- 
conserving diffusion Monte Carlo [2], In order to reproduce accurately the equation of state fully convergent 
ab initio calculations required for the reference potential surface. In Mark I model we fit to empirical 
potential [1], which effectively included manybody interactions, and as a result, the experimental density 
was reproduced when dispersion interactions were scaled up by 37%. 



FIG. 1: QDO model for water. 

The system is thus defined by monomers with “perfect” (by construction) long-range electrostatics and 
all resulting condensed phase behavior emerges as a model prediction. No additional thermodynamic data 
or liquid state information was used. 

The coarse grained Hamiltonian for the quantum Drude oscillator (QDO) water model is 

h = J2 T i Iis) + </> (Coul) (R) + E ^ (rep) CHoij)+ 4 d) (r), 

i j+i 

H( D ^o(r,R) = ^ D) (R)^ 0 (r,R), 

H (D) = E (D - Rci) 2 ) + </> (Coul) (r, R). 

Here, T-'" gi is the classical rigid body kinetic energy of water molecule i; 9A r -vcctor R represents the co¬ 
ordinates of all the water molecules in the system; (/^ c:oul) (R) is the intermolecular Coulomb interaction 
energy between the fixed monomer gas phase charge distributions only; </>( rep ) (Roij) is an oxygen-centered 
pair-wise repulsion. The quantities E^ D> (R) and -00 ( r : R) are the ground state Born-Oppenheimer energy 
surface and wavefunction of the QDO electronic structure, respectively; r j and T, are the position and quan¬ 
tum kinetic energy of the oscillator centered on molecule i. The oscillators centered at R c , are characterized 
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by parameters {p,u,q}. Finally, 0( Coul * fr. R) is intermolecular Coulomb interaction between the Drude 
particles of charge q~, at position r, their center of oscillation with charge q + , at position R c , and the fixed 
monomer charge distribution. Regularization of the Coulomb interaction is discussed elsewhere [3]. 

Finite temperature condensed phase calculations were performed using adiabatic path integral molecular 
dynamics for quantum Drude oscillators (APIMD-QDO) [4-6] in the canonical and isothermal-isobaric 
ensembles [7]. We used the discretized (high-temperature) path integral approximation to the partition 
function with the number of beads (the Trotter number) P = 96 and adiabaticity factor 7 = 16 and close to 
Widom line densities, 7 = 32. 

All calculations were performed for canonical ( NVT ) or isothermal-isobaric ( NpT ) conditions using 
cubic systems containing 300 water molecules with 3D periodic boundary conditions and Ewald summa¬ 
tion technique for electrostatic interactions. At 673 K we recalculated all density points with N = 1000 
molecules and the results were statistically identical to that smaller system. This is the consequence of the 
absence of any truncation scheme in intermolecular interactions and of the short range of structural and 
orientational ordering at high temperatures. The integration timestep for adiabaticity factor 7 = 16 was 
dt = 0.15 fs (6.25 Ti/Eh) and a typical integration time was 150 ps (lif’dt), which requires ca. 27 h of CPU 
time on BlueGene/Q running on 8 K cores. 

Finite size effects estimation 

In order to verify that our results converged with respect to integration time and system size, for one 
thermodynamic point, T* = 1.036 and density p* = 0.82 we performed a set of test calculations using 
SWM4-NDP model for systems with N = 300, 1,000, 4,000, 12,500, and 100,000 molecules. The Drude 
oscillator degrees of freedom were treated as dynamic variabled and extended Lagrangian equations of 
motion were integrated with 0.5 fs timestep, and the temperature was controlled by two thermostats with 
classical Drude oscillators kept close to 1 K by Langevin thermostat. The obtained results for the dipole 
moment convergence are shown in Fig. 2. They fit well to 1/N (shown in the Figure by the dashed line) 
and from this fit the dipole moment extrapolated to infinite size is only than 0 . 1 % higher than the value 
obtained with the smallest studied system. No systematic drift in statistical properties has been observed 
after equilibrating the systems for at least 200 ps. 

Special care was taken to ensure that there are no systematic drift in computer properties, for which 
long runs, from 0.1 ps for N = 300 to 1 ns for N = 10 5 have been performed. In all cases pressure was 
p = 24.81(4) MPa. 
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FIG. 2: System size dependence of the molecular dipole moment calculated for the classical Drude model at reduced 
density p* = 0.82 and reduced temperature T* = 1.036. The error bars correspond to two standard deviations (95% 
confidence level). 


Hydrogen bonding 

When calculating various HB distributions we used the geometric R — (3 distance-angle definition of the 
hydrogen bond [8], in which two water molecules are considered hydrogen-bonded if and only if the 0-0 
distance between two molecules is shorter than the certain threshold value, i?HB, and the angle (3 between 
O-H bond and the candidate molecule oxygen, O — H ■ ■ ■ O, is smaller than the specified criterion /3hb- We 
have found that the values Ehb = 3.5 A and /3hb = 30° are suitable for the studied cases (see Fig. 3). 


Molecular moment scaling 


Molecular dipole moment scales linearly with density. This could be understood from simple argument. 
The second contribution can be treated in mean-field approximation. By the chain rule for differentia- 


d\p\ = d\p\dE^dR 

dp 3Ei dR dp' K ’ 

Here, E) is the projection of the local Coulomb field on the molecule along the ZHOH bisector, and R 
is a measure of average molecular separation. The first factor on the RHS of (1) is the molecular po¬ 
larizability a along the field direction, to leading order. Defining the density as p = rail " 3 we have 
dR/dp ~ —/i > 4 /m. For the local field dependence on R, we assume leading term should be due to the 
permanent dipole po, the separation R and an orientational factor f(0) according to pnf(O ') /i? 3 . Therefore 
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FIG. 3: Potential of mean force (in kT units) map in 0-0 and ZOH ■ ■ • O dimensions at various densities along the 
673 K isotherm. The dashed lines of contrast colors fence the area of hydrogen-bonded molecules using our criterion. 


( dE/dR ) ~ — fj,of(0)/R 4 and, finally, dfi/dp ~ a p,of(9)/m. Consequently, we expect the molecular 
dipole to depend approximately linearly on density with the observed change of slope indicative of an in¬ 
crease in local density as the HB network coalesces. This is consistent with the density dependence of the 
Onsager reaction field [9]. 
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FIG. 4: Dipole moment correlation function for the first solvation shell for the Stockmayer fluid, left column; QDO 
water, right column at T* = 1,0367' c at denoted densities. 
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Polarizable Stockmayer model 

The original Stockmayer model, consisting of the Lennard-Jones (LJ) isotropic part and embedded 
point dipole, provides one of the simplest models for polar liquids. In particular, we study a dipole po¬ 
larizable version using the parameters of Ref. [10] as a baseline. In order to run the simulations in the 
NAMD code [11], we replaced the point dipole by a pah - of ±0.52|e| point charges 0.6 A apart, giving 
the permanent dipole moment of 1.5 D. A classical Drude model embedded to generate induction in the 
dipole limit. The oscillator consists of a point charge of —0.5|e| harmonically tethered to the midpoint 
between the charges forming the fixed dipole; the latter bears a neutralizing charge of +0.5|e|. The gener¬ 
ally accepted value for the water polarizability, a = 1.444 A 1 2 3 4 5 6 7 8 9 10 11 , defines the classical Drude spring constant, 
k = 240 kJ mol -1 A 2 . Apart from the Lennard-Jones terms between the Drude tether sites, with parame¬ 
ters e = 3.1785 kJ/mol and a = 2.955 A, a regularising LJ repulsion is added between the oscillators with 
parameters e = 0.36854 kJ/mol, a = 2.955 A. 

Association topology 

In order to study the short-range orientational correlations of water and the Stockmayer fluid, we calcu¬ 
lated the probability distrinution function for two angles 9i i E {1,2} between molecular dipole moments 
and the vector connecting their - COMs for the molecules in the first solvation shell (as defined by the 1st 
minima in goo(r) radial distribution functions). The results presented in Fig. 4 show that while water 
prefers tetrahedral coordination, the Stockmayer fluid adopts a linear association topology. 
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